{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adaptivity\n",
    "==="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the geometry by 2D Netgen-OpenCascade modeling (new in Netgen/NGSolve 2105):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "\n",
    "def MakeGeometryOCC():\n",
    "    base = Rectangle(1, 0.6).Face()\n",
    "    chip = MoveTo(0.5,0.15).Line(0.15,0.15).Line(-0.15,0.15).Line(-0.15,-0.15).Close().Face()\n",
    "    top = MoveTo(0.2,0.6).Rectangle(0.6,0.2).Face()\n",
    "    base -= chip\n",
    "\n",
    "    base.faces.name=\"base\"\n",
    "    chip.faces.name=\"chip\"\n",
    "    chip.faces.col=(1,0,0)\n",
    "    top.faces.name=\"top\"\n",
    "    geo = Glue([base,chip,top])\n",
    "    geo.edges.name=\"default\"\n",
    "    geo.edges.Min(Y).name=\"bot\"\n",
    "    DrawGeo(geo)\n",
    "    return OCCGeometry(geo, dim=2)\n",
    "\n",
    "geo = MakeGeometryOCC()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the geometry by curves (old-style segments geometry):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   point numbers 0, 1, ... 11\n",
    "#   sub-domain numbers (1), (2), (3)\n",
    "#  \n",
    "#\n",
    "#             7-------------6\n",
    "#             |             |\n",
    "#             |     (2)     |\n",
    "#             |             |\n",
    "#      3------4-------------5------2\n",
    "#      |                           |\n",
    "#      |             11            |\n",
    "#      |           /   \\           |\n",
    "#      |         10 (3) 9          |\n",
    "#      |           \\   /     (1)   |\n",
    "#      |             8             |\n",
    "#      |                           |\n",
    "#      0---------------------------1\n",
    "#\n",
    "\n",
    "def MakeGeometry():\n",
    "    from netgen.geom2d import SplineGeometry\n",
    "    geometry = SplineGeometry()\n",
    "    \n",
    "    # point coordinates ...\n",
    "    pnts = [ (0,0), (1,0), (1,0.6), (0,0.6), \\\n",
    "             (0.2,0.6), (0.8,0.6), (0.8,0.8), (0.2,0.8), \\\n",
    "             (0.5,0.15), (0.65,0.3), (0.5,0.45), (0.35,0.3) ]\n",
    "    pnums = [geometry.AppendPoint(*p) for p in pnts]\n",
    "    \n",
    "    # start-point, end-point, boundary-condition, left-domain, right-domain:\n",
    "    lines = [ (0,1,\"bot\",1,0), (1,2,\"outer\",1,0), (2,5,\"outer\",1,0), (5,4,\"inner\",1,2), (4,3,\"outer\",1,0), (3,0,\"outer\",1,0), \\\n",
    "              (5,6,\"outer\",2,0), (6,7,\"outer\",2,0), (7,4,\"outer\",2,0), \\\n",
    "              (8,9,\"inner\",3,1), (9,10,\"inner\",3,1), (10,11,\"inner\",3,1), (11,8,\"inner\",3,1) ]\n",
    "        \n",
    "    for p1,p2,bc,left,right in lines:\n",
    "        geometry.Append([\"line\", pnums[p1], pnums[p2]], bc=bc, leftdomain=left, rightdomain=right)\n",
    "\n",
    "    geometry.SetMaterial(1,\"base\")\n",
    "    geometry.SetMaterial(2,\"top\")    \n",
    "    geometry.SetMaterial(3,\"chip\")\n",
    "\n",
    "    return geometry\n",
    "\n",
    "# geo = MakeGeometry()\n",
    "# Draw(geo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Piece-wise constant coefficients in sub-domains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mesh = Mesh(geo.GenerateMesh(maxh=0.2))\n",
    "\n",
    "fes = H1(mesh, order=3, dirichlet=\"bot\", autoupdate=True)\n",
    "u, v = fes.TnT()\n",
    "\n",
    "lam = CoefficientFunction([1, 1000, 10])\n",
    "a = BilinearForm(fes)\n",
    "a += lam*grad(u)*grad(v)*dx\n",
    "\n",
    "# heat-source in inner subdomain\n",
    "f = LinearForm(fes)\n",
    "f += 1*v*dx(definedon=\"chip\")\n",
    "\n",
    "c = Preconditioner(a, type=\"multigrid\", inverse=\"sparsecholesky\")\n",
    "\n",
    "gfu = GridFunction(fes, autoupdate=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble and solve problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SolveBVP():\n",
    "    a.Assemble()\n",
    "    f.Assemble()\n",
    "    inv = CGSolver(a.mat, c.mat)\n",
    "    gfu.vec.data = inv * f.vec\n",
    "    \n",
    "SolveBVP()\n",
    "Draw (gfu, mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gradient recovery error estimator: Interpolate finite element flux \n",
    "\n",
    "$$\n",
    "q_h := I_h (\\lambda \\nabla u_h)\n",
    "$$\n",
    "\n",
    "and take difference as element error indicator:\n",
    "\n",
    "$$\n",
    "\\eta_T := \\tfrac{1}{\\lambda} \\| q_h - \\lambda \\nabla u_h \\|_{L_2(T)}^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = []    # l = list of estimated total error\n",
    "space_flux = HDiv(mesh, order=2, autoupdate=True)\n",
    "gf_flux = GridFunction(space_flux, \"flux\", autoupdate=True)\n",
    "\n",
    "def CalcError():\n",
    "    \n",
    "    # FEM-flux \n",
    "    flux = lam * grad(gfu)        \n",
    "    # interpolate into H(div)\n",
    "    gf_flux.Set(flux) \n",
    "    \n",
    "    # compute estimator:\n",
    "    err = 1/lam*(flux-gf_flux)*(flux-gf_flux)\n",
    "    eta2 = Integrate(err, mesh, VOL, element_wise=True)\n",
    "    l.append ((fes.ndof, sqrt(sum(eta2))))\n",
    "    print(\"ndof =\", fes.ndof, \" toterr =\", sqrt(sum(eta2)))\n",
    "    \n",
    "    # mark for refinement:\n",
    "    maxerr = max(eta2)\n",
    "    # marking with Python loop:\n",
    "    # for el in mesh.Elements():\n",
    "    #    mesh.SetRefinementFlag(el, eta2[el.nr] > 0.25*maxerr)\n",
    "    \n",
    "    # marking using numpy:\n",
    "    mesh.ngmesh.Elements2D().NumPy()[\"refine\"] = \\\n",
    "        eta2.NumPy() > 0.25*maxerr\n",
    "  \n",
    "CalcError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adaptive loop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "level = 0\n",
    "while fes.ndof < 50000:  \n",
    "    mesh.Refine()\n",
    "    SolveBVP()\n",
    "    CalcError()\n",
    "    level = level+1\n",
    "    if level%5 == 0:\n",
    "        Draw (gfu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Draw (gfu);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(\"ndof\")\n",
    "plt.ylabel(\"H1 error-estimate\")\n",
    "ndof,err = zip(*l)\n",
    "plt.plot(ndof,err, \"-*\")\n",
    "\n",
    "plt.ion()\n",
    "plt.show();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
